%%%%%%%%%%%%%%%%%%
% options market making model, slight adjustment of return factor model
% Bas = Structural factors (HC, VR, JR)+[HRP, VRP] + Option option interest 
% Liuren Wu, liuren.wu@baruch.cuny.edu
%%%%%%%%%%%%%%%%%%
function [Af,Bf,R2s,stats,VB] =FunDailyORSestimationBAS(t,y,x,bas,HVV,IV,ids,dds,dd,nlag,nx,winth,minN)

tt=dds==dd(t);  sdt=ids(tt);
yt=y(tt);
bast=winsorize(bas(tt),winth);

Nt=sum(tt);
rsv=NaN(Nt,nlag);
jj=1:min(t-1,nlag);
nj=length(jj);
Biv=NaN(3,nj);Bi0=ones(3,1)/3;
for j=1:nj
    ll=jj(j);
    ttp=dds==dd(t-ll);sdp=ids(ttp);
    ytp=y(ttp);
    [~,ind1,ind2]=intersect(sdt,sdp);
    rsv(ind1,j)=ytp(ind2);
    
    HVVp=HVV(ttp,:);
    
    tti=dds==dd(t-ll+1);HVi=HVV(tti,1);sdi=ids(tti);
    [~,ind1,ind2]=intersect(sdi,sdp);
    yi=HVi(ind1);
    Xi=HVVp(ind2,:);
    ii=isfinite(yi+sum(Xi,2)); yi=yi(ii);Xi=Xi(ii,:);
    XX=Xi'*Xi;
    if rank(XX)<3
        Pi=diag(diag(XX))*1e-4;
        Biv(:,j)=(XX+Pi)\(Xi'*yi + Pi*Bi0);
    else
        Biv(:,j)=Xi\yi;
    end
end
Bi=nanmean(Biv,2);

HVVt=HVV(tt,:);FHV=HVVt*Bi;
VRP=100*(IV(tt) -FHV); % VRP
HRP=nanmean(rsv(:,2:nlag),2);%HRP

tvrx=(nanstd(rsv,0,2)); %rolling-window return volatility estimator

xt=[x(tt,:),HRP,VRP]; 


xt=winsorize(xt,winth);
tvr=winsorize(tvrx,winth);
%1. Descriptor stats
Zt=[bast,xt,tvr];
stats=[...
    nanmean(Zt);
    nanstd(Zt,0);
    prctile(Zt,[10,25,50,75,90]');
    skewness(Zt,0,1);
    kurtosis(Zt,0,1)-3;
    100*corr(yt,Zt,'rows','pairwise')];

% %regression
SZt=NaN(Nt,nx); 
indfx=zeros(nx,1);
SXV=NaN(1,nx); %scale
SZ2=SZt;
for j=1:nx
    indf=isfinite(xt(:,j)+yt); Nj=sum(indf);
    if Nj>minN
        indfx(j)=1;
        xtj=xt(indf,j);sxj=std(xtj);mxj=mean(xtj);
        if j<=3
            SZt(indf,j)=xtj./sxj;
        else
            SZt(indf,j)=(xtj-mxj)./sxj;
        end
        SXV(j)=sxj;
        SZ2(indf,j)=xtj;
    end
end

Bf=NaN(nx,1); R2s=NaN(1,1);Af=NaN(1,1); VB=Bf;
indff=find(indfx==1);
if ~isempty(indff)
    Ft = SZt(:,indff);     [Af(1),Bf(indff,1),R2s,VB]=FunFMRegression(Ft,bast,SZ2);
end

end


function [A,B,AR2,vb,CO]=FunFMRegression(Ft,yt,SZ2)
indf=isfinite(sum(Ft,2)+yt); Nf=sum(indf);

Xtf=[ones(Nf,1),Ft(indf,:)];ytf=yt(indf); 
B=Xtf\ytf; 
yh=Xtf*B; 
ef=ytf-yh;
nx=size(Xtf,2);
AR2= 1-(sum(ef.^2)./(Nf-nx))/var(ytf);

A=B(1); B=B(2:end);
Stf=SZ2(indf,:);
C=Xtf\Stf;
CO=diag(C(2:end,:));

vb=B.*(cov(Xtf(:,2:end))*B)./var(ytf);



end




